################################################################################################

options(stringsAsFactors=F)
library(parallel)
library(data.table)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--cluster"), type = "character"),
    make_option(c("--batch"), type = "character"),
    make_option(c("--barcode"), type = "character"),
    make_option(c("--peak"), type = "character"),
    make_option(c("--frag"), type = "character"),
    make_option(c("--output"), type = "character")
)

###########################################################################################

if(1!=1){

	clustern <- 15
	batchn <- 2
	barcode_dir <- "/public/home/xxf2019/20231121_singleMuti/results/cluster_all_result/cell_fragment"
	peak_file <- "/public/home/xxf2019/20231121_singleMuti/results/cluster_all_result/cluster15/consensus/cluster15.fwp.filter.non_overlapping.bed"
	frag_dir <- "/public/home/xxf2019/20231121_singleMuti/results/cluster_all_result/cluster15/fragment"
	output_dir <- "/public/home/xxf2019/20231121_singleMuti/results/cluster_all_result/cluster15/count"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

clustern <- opt$cluster
batchn <- opt$batch
barcode_dir <- opt$barcode
peak_file <- opt$peak
frag_dir <- opt$frag
output_dir <- opt$output

################################################################################################

testis_cell <- read.csv(paste0(barcode_dir,"/testis0",batchn,"_all_qc_barcode.csv"))

## function ##
count_func <- function(cell,peakt,frag){
	countm <- data.frame(rbindlist(mclapply(1:nrow(cell),function(celli){
		print(celli)
		in_count <- subset(frag,V10==cell$V1[celli])[,c("V4","V11")]

		if(nrow(in_count)!=0){
			colnames(in_count) <- c("peak","count")
			in_count_cell <- data.frame(table(in_count$peak))
			in_count_mc <- subset(in_count_cell,Freq>1)
			if(nrow(in_count_mc)>0){
				#print(1)
				frag_count1 <- data.frame(rbindlist(lapply(1:nrow(in_count_mc),function(peaki){
					peakp <- subset(in_count,peak==in_count_mc$Var1[peaki])
					return(data.frame(peak=in_count_mc$Var1[peaki],count=sum(as.numeric(peakp$count))))
				})))
				frag_count2 <- subset(in_count,peak%in%subset(in_count_cell,Freq==1)$Var1)
				frag_count3 <- data.frame(peak=subset(peakt,!V4%in%in_count_cell$Var1)$V4,count=0)
				frag_count <- rbind(frag_count1,frag_count2,frag_count3)
				cell_bid <- cell$cellID[celli]
				colnames(frag_count) <- c("peak",cell_bid)
				frag_count <- data.frame(t(frag_count))
				colnames(frag_count) <- frag_count[1,]

			}else{
				#print(2)
				frag_count2 <- subset(in_count,peak%in%subset(in_count_cell,Freq==1)$Var1)
				frag_count3 <- data.frame(peak=subset(peakt,!V4%in%in_count_cell$Var1)$V4,count=0)
				frag_count <- rbind(frag_count2,frag_count3)
				cell_bid <- cell$cellID[celli]
				colnames(frag_count) <- c("peak",cell_bid)
				frag_count <- data.frame(t(frag_count))
				colnames(frag_count) <- frag_count[1,]
			}
		}else{
			#print(3)
			frag_count <- data.frame(peak=subset(peakt,!V4%in%in_count$Var1)$V4,count=0)
			cell_bid <- cell$cellID[celli]
			colnames(frag_count) <- c("peak",cell_bid)
			frag_count <- data.frame(t(frag_count))
			colnames(frag_count) <- frag_count[1,]
		}
		frag_count <- frag_count[,peakt$V4]
		return(frag_count)

	},mc.cores=20)))

	countm
}

################################################################################################
# count #
cluster_peak <- read.table(peak_file)

testis_frag <- data.frame(fread(paste0(frag_dir,"/cluster",clustern,"_testis0",batchn,"_frag_overlap.txt"),header=FALSE))
cluster_testis_count <- count_func(testis_cell,cluster_peak,testis_frag)
cluster_testis_count <- cluster_testis_count[-seq(1,nrow(cluster_testis_count),2),]
cluster_testis_count <- data.frame(t(cluster_testis_count))
colnames(cluster_testis_count) <- testis_cell$cellID
save(cluster_testis_count,file=paste0(output_dir,"/cluster",clustern,"_testis0",batchn,"_count.Rdata"))


cell=testis_cell
peakt=cluster_peak
frag=testis_frag